
#include "geometry.h"
#include "utility.hpp"
#include <fstream>
#include <iomanip>
#include <iostream>
using namespace std;

Geometry::Geometry(const Config &usrInput)
{
  ReadBoundFile(usrInput.GetBoundFname());

  // 可以用 map 和多态
  switch (usrInput.GetMeshFormat())
  {
  case MeshFormat_t::ASSCII:
    ReadMeshASSCII(usrInput.GetMeshFname());
    break;
  case MeshFormat_t::Binary:

    break;
  default:
    cerr << "\nmesh format err\n"
         << endl;
    // causes normal program termination with cleaning up
    exit(EXIT_FAILURE);
    break;
  }
  // 配置边界信息
  SetBoundDescrpt();
  cutExchgVol.SetData(i2, j2);
  // InitMetrics();
}

void Geometry::InitMetrics(void)
{
  // 计算了所有点构成的面 包含部分虚单元的面
  for (int j = 0; j <= jl; ++j)
    for (int i = 0; i <= imax; ++i)
    {
      si(i, j).sx = coord(i, j).y - coord(i, j + 1).y;
      si(i, j).sy = coord(i, j + 1).x - coord(i, j).x;
    }
  for (int j = 0; j <= jmax; ++j)
    for (int i = 0; i <= il; ++i)
    {
      sj(i, j).sx = coord(i + 1, j).y - coord(i, j).y;
      sj(i, j).sy = coord(i, j).x - coord(i + 1, j).x;
    }

  for (int j = 2; j <= j2; ++j)
    for (int i = 2; i <= i2; ++i)
    {
      // 当 i j 按右手螺旋  拇指指向屏幕外 体积为正 若为屏幕内，为负
      vol(i, j) = fabs(0.5 * ((coord(i, j).x - coord(i + 1, j + 1).x) *
                                  (coord(i + 1, j).y - coord(i, j + 1).y) +
                              (coord(i, j + 1).x - coord(i + 1, j).x) *
                                  (coord(i, j).y - coord(i + 1, j + 1).y)));
      // cout << vol(i, j) << "\n";
    }

  // ! ============================================================================
  vol(0, 0) = 0.0;
  vol(0, 1) = 0.0;
  vol(1, 0) = 0.0;
  vol(1, 1) = 0.0;
  vol(il, 0) = 0.0;
  vol(il, 1) = 0.0;
  vol(imax, 0) = 0.0;
  vol(imax, 1) = 0.0;
  vol(0, jl) = 0.0;
  vol(0, jmax) = 0.0;
  vol(1, jl) = 0.0;
  vol(1, jmax) = 0.0;
  vol(il, jl) = 0.0;
  vol(il, jmax) = 0.0;
  vol(imax, jl) = 0.0;
  vol(imax, jmax) = 0.0;
  // ! ============================================================================
  // 对虚单元体积赋值
  for (int i = 2; i <= i2; ++i)
  {
    vol(i, 0) = vol(i, 1) = vol(i, 2);
    vol(i, jmax) = vol(i, jl) = vol(i, j2);
  }

  for (int j = 2; j <= j2; ++j)
  {
    vol(0, j) = vol(1, j) = vol(2, j);
    vol(imax, j) = vol(il, j) = vol(i2, j);
  }

  //================================
  CutExchg<real_t> boundCut;
  for (int i = 0; i < nBound; ++i)
  {
    if (BoundDescrp[i].BType == Bound_t::InnerBound)
    // (BoundDescrp[i].BType == Bound_t::Periodic))
    {
      cutExchgVol(BoundDescrp[i], vol);
    }
  }
  //================================

  // 检查体积和面积
  real_t smax{-1.0e10}, volMin{1.0e10}, volMax{-1.0e10};
  real_t sumx, sumy, sumNorm;
  for (int j = 2; j <= j2; ++j)
    for (int i = 2; i <= i2; ++i)
    {
      sumx = si(i, j).sx - si(i + 1, j).sx + sj(i, j).sx - sj(i, j + 1).sx;
      sumy = si(i, j).sy - si(i + 1, j).sy + sj(i, j).sy - sj(i, j + 1).sy;
      sumNorm = sqrt(sumx * sumx + sumy * sumy);
      smax = fmax(sumNorm, smax);
      volMin = fmin(volMin, vol(i, j));
      volMax = fmax(volMax, vol(i, j));
    }
  cout << scientific << setprecision(15);
  cout << "\n max sum(area) = " << smax << "\n min(vol) = " << volMin
       << "\n max(vol) = " << volMax << "\n";
  // cout << fmax(2.0,1.0)  << endl;
}

// 考虑的只是单块  多块需要大改
void Geometry::ReadBoundFile(const std::string &boundFname)
{
  cout << "start reading boundary file\n";

  fstream boundF{boundFname};
  if (!boundF)
  {
    cerr << "open boundary file err\n"
         << endl;
    exit(EXIT_FAILURE);
  }
  int tmp, blockNum;
  string faceName;

  boundF >> tmp;
  boundF >> blockNum;
  if (blockNum != 1)
  {
    cerr << "only support one block now \n"
         << endl;
    exit(EXIT_FAILURE);
  }

  boundF >> nci >> ncj;
  // 输入的是网格点数
  --nci;
  --ncj;
  boundF >> faceName;
  boundF >> nBound;
  BoundDescrp.resize(nBound);

  // 不是想封装，而是想降维， 需要使用vector 动态添加
  struct TmpBound
  {
    int ib, ie, jb, je;
    Bound_t boundType;
  };

  vector<TmpBound> genericBound;
  TmpBound tmpB;
  int btype;
  // 读取 generic 边界文件格式
  for (int i = 0; i < nBound; ++i)
  {
    boundF >> tmpB.ib >> tmpB.ie >> tmpB.jb >> tmpB.je >> btype;
    tmpB.boundType = static_cast<Bound_t>(btype);
    // 不会出现负的索引  负的索引只是用于配对  二维不需要配对
    // 网格点与网格单元的关系转换
    // 0 1 是虚单元 修正起始索引
    ++tmpB.ib;
    ++tmpB.ie;
    ++tmpB.jb;
    ++tmpB.je;

    genericBound.push_back(tmpB);

    if (tmpB.boundType == Bound_t::InnerBound)
      --i;
  }

  // 将 generic 边界格式转换为 blazik 中的格式
  int count{0};
  for (int i = 0; i < genericBound.size(); ++i)
  {
    tmpB = genericBound[i];
    // 物理边界
    if (tmpB.boundType != Bound_t::InnerBound)
    {
      BoundDescrp_t &mid = BoundDescrp[count];
      mid.BType = tmpB.boundType;

      // i 面
      if (tmpB.ib == tmpB.ie)
      {
        mid.Bn = (tmpB.ib == 2) ? 1 : 3;
        mid.beg = tmpB.jb;
        mid.end = tmpB.je;
        // 保证了物理边界索引从小到大
        if (mid.end < mid.beg)
          swap(mid.beg, mid.end);
        // 修正为单元索引
        --mid.end;

        mid.Bns = -10;
        mid.begs = mid.ends = -10;
      }
      else
      // j 面
      {
        mid.Bn = (tmpB.jb == 2) ? 2 : 4;
        mid.beg = tmpB.ib;
        mid.end = tmpB.ie;
        // 保证了物理边界索引从小到大
        if (mid.end < mid.beg)
          swap(mid.beg, mid.end);
        --mid.end;

        mid.Bns = -10;
        mid.begs = mid.ends = -10;
      }

      ++count;
      // 从0起始
      --mid.Bn;
    }
    else
    // ==================================================
    // 内边界  内边界仍然需要判断 beg 是否 大于 end  谨慎起见
    {
      BoundDescrp_t &mid = BoundDescrp[count];
      mid.BType = tmpB.boundType;
      // 当前面为 i 面
      if (tmpB.ib == tmpB.ie)
      {
        mid.Bn = (tmpB.ib == 2) ? 1 : 3;
        mid.beg = tmpB.jb;
        mid.end = tmpB.je;
        // 修正为单元索引
        if (mid.end < mid.beg)
          --mid.beg;
        else
          --mid.end;
      }
      else
      // 当前面为 j 面
      {
        mid.Bn = (tmpB.jb == 2) ? 2 : 4;
        mid.beg = tmpB.ib;
        mid.end = tmpB.ie;
        // 修正为单元索引
        if (mid.end < mid.beg)
          --mid.beg;
        else
          --mid.end;
      }
      // 邻接面 “源面”
      ++i;
      tmpB = genericBound[i];
      // 紧接着下面的一行最后一个为块号
      // 邻接块为 i 面
      if (tmpB.ib == tmpB.ie)
      {
        mid.Bns = (tmpB.ib == 2) ? 1 : 3;
        mid.begs = tmpB.jb;
        mid.ends = tmpB.je;

        if (mid.ends < mid.begs)
          --mid.begs;
        else
          --mid.ends;
      }
      else
      // 邻接块为 j 面
      {
        mid.Bns = (tmpB.jb == 2) ? 2 : 4;
        mid.begs = tmpB.ib;
        mid.ends = tmpB.ie;

        if (mid.ends < mid.begs)
          --mid.begs;
        else
          --mid.ends;
      }
      ++count;
      --mid.Bn;
      --mid.Bns;
    }
  }

  if (count != nBound)
  {
    cerr << "generic boundary format to blazik format err\n"
         << endl;
    exit(EXIT_FAILURE);
  }
  boundF.close();
  cout << "read boundary file success" << endl;
}

void Geometry::ReadMeshASSCII(const std::string &meshFname)
{
  cout << "start reading mesh file\n";

  ifstream meshF{meshFname};
  if (!meshF)
  {
    cerr << "open mesh file err\n"
         << endl;
    exit(EXIT_FAILURE);
  }

  // plot3d 网格文件要以多块输出
  int blockNum, ni, nj, nk;
  meshF >> blockNum;
  if (blockNum != 1)
  {
    cerr << "only support one block now\n"
         << endl;
    exit(EXIT_FAILURE);
  }
  meshF >> ni >> nj >> nk;
  if ((--ni != nci) && (--nj != ncj) && (nk != 1))
  {
    cerr << "mesh file content err\n"
         << endl;
    exit(EXIT_FAILURE);
  }

  imax = nci + 3;
  jmax = ncj + 3;
  il = imax - 1;
  jl = jmax - 1;
  i2 = imax - 2;
  j2 = jmax - 2;
  // 跟 blazik 同步 后面有需求再改
  coord.ReConstruct(imax + 1, jmax + 1);
  si.ReConstruct(imax + 1, jmax + 1);
  sj.ReConstruct(imax + 1, jmax + 1);
  vol.ReConstruct(imax + 1, jmax + 1);

  for (int j = 2; j <= jl; ++j)
    for (int i = 2; i <= il; ++i)
    {
      meshF >> coord(i, j).x;
    }

  for (int j = 2; j <= jl; ++j)
    for (int i = 2; i <= il; ++i)
    {
      meshF >> coord(i, j).y;
    }
// 为了输出查看网格
#ifdef TEST
  for (int j = 2; j <= jl; ++j)
    for (int i = 2; i <= il; ++i)
    {
      meshF >> coord(i, j).z;
    }
#endif
  meshF.close();

  // 虚单元坐标赋值
  for (int i = 2; i < imax; ++i)
  {
    coord(i, 0) = coord(i, 1) = coord(i, 2);
    coord(i, jmax) = coord(i, jl);
  }
  for (int j = 0; j <= jmax; ++j)
  {
    coord(0, j) = coord(1, j) = coord(2, j);
    coord(imax, j) = coord(il, j);
  }

  cout << "read mesh file success" << endl;
}

void Geometry::SetBoundDescrpt(void)
{
  const int dum[4][2]{{1, 0}, {1, 0}, {il, imax}, {jl, jmax}};
  const int ins[4][2]{{2, 3}, {2, 3}, {i2, i2 - 1}, {j2, j2 - 1}};
  // 后面加的
  for (int i = 0; i < nBound; ++i)
  {
    BoundDescrp_t &tmp = BoundDescrp[i];
    tmp.dum1 = dum[tmp.Bn][0];
    tmp.dum2 = dum[tmp.Bn][1];
    tmp.ins1 = ins[tmp.Bn][0];
    tmp.ins2 = ins[tmp.Bn][1];
    // 保证当前单元边界索引从小到大
    if (tmp.beg > tmp.end)
    {
      swap(tmp.beg, tmp.end);
      swap(tmp.begs, tmp.ends);
    }
    tmp.isIface = ((tmp.Bn == 0) || (tmp.Bn == 2)) ? true : false;
    tmp.isIfaceAdj = ((tmp.Bns == 0) || (tmp.Bns == 2)) ? true : false;
    tmp.lsteps = (tmp.begs < tmp.ends) ? 1 : -1;
  }
}

#ifdef TEST
void Geometry::PrintMeshASSCII(const std::string &outMeshFname)
{
  ofstream meshF{outMeshFname};
  meshF << 1 << "\n";
  meshF << nci + 1 << "\n"
        << ncj + 1 << "\n"
        << 1 << "\n";
  meshF << scientific << setprecision(15);
  for (int k = 0; k < 3; ++k)
  {
    for (int j = 2; j <= jl; ++j)
    {
      for (int i = 2; i <= il; ++i)
      {
        switch (k)
        {
        case 0:
          meshF << coord(i, j).x << "\t";
          break;
        case 1:
          meshF << coord(i, j).y << "\t";
          break;
        case 2:
          meshF << coord(i, j).z << "\t";
          break;
        }
        if (i % 4 == 0)
          meshF << "\n";
      }
    }
    meshF << "\n";
  }
  meshF.close();
}
#endif

void Geometry::PrintBoundary(ostream &os) const
{
  os << "num cells i = " << nci << "\n";
  os << "num cells j = " << ncj << "\n\n";

  os << "nBound = " << nBound << "\n";
  for (int i = 0; i < nBound; ++i)
  {
    const BoundDescrp_t &tmp = BoundDescrp[i];
    os << static_cast<int>(tmp.BType) << "\t" << tmp.Bn << "\t" << tmp.beg
       << "\t" << tmp.end << "\t" << tmp.Bns << "\t" << tmp.begs << "\t"
       << tmp.ends << "\t"
       << "   dum2 = " << tmp.dum2 << "   dum1 = " << tmp.dum1
       << "   ins1 = " << tmp.ins1 << "   ins2 = " << tmp.ins2 << "\n";
  }
}